Monte Carlo discretization of general relativistic radiation transport 
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An indirect, hybrid Monte Carlo discretization of general relativistic kinetic theory suitable for the 
development of numerical schemes for radiation transport is presented. The discretization is based 
on surface flux estimators obtained from a local decomposition of the distribution function, and can 
handle optically thick regions by means of formal solutions within each cell. Furthermore, the scheme 
is designed for parallel implementation, and it admits the use of adaptive techniques by virtue of 
leaving all probability density functions unspecified. Some considerations for numerical uses of the 
scheme are discussed. 
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I. INTRODUCTION 

Radiation transport is an important process in many 
astrophysical systems. In those models which involve 
compact objects, a notion which encompasses active 
galactic nuclei, supernovae, collapsars, X-ray binaries 
and microquasars, and compact binary mergers as mod- 
els for short-hard gamma-ray bursts, either photons or 
neutrinos can have a role as a dynamical agent, influ- 
ence the chemical and nuclear composition of the mat- 
ter, and give rise to directly observable signals. 

In recent years, numerical techniques for handling 
general relativistic models, in particular including those 
involving black holes and relativistic magnetized fluids, 
have made substantial progress H H HM S IE H H U 
Hot [Til [l2L [l3l [Till . Current numerical codes are able 
to evolve isolated, accreting and binary black holes and 
neutron stars for long times compared to the dynamical 
time scales of the compact objects. 

Most models, however, do not account for radiation 
transport. In some cases this is a reasonable approxima- 
tion, particularly when the radiation is trapped and only 
influences the system on diffusion timescales, but there 
are many interesting situations where the interaction of 
matter with radiation fields cannot be neglected. Radi- 
ation can remove or supply energy and momentum to 
the matter, and thus influence the local thermodynam- 
ical state and composition of the material. In turn, this 
affects the bulk dynamical properties of fluid elements 
in a numerical simulation. 

A major difficulty with discretizing radiation trans- 
port is its high-dimensional nature. Each radiation field 
depends on three positional and three momentum vari- 
ables, which makes even storing the full field state at 
a particular time, and with enough resolution, a chal- 
lenging proposition. One popular strategy is to make 
the diffusion assumption 111511 with additional flux limiters 
to handle the regions near free-streaming fl^H . Another 
strategy reduces the problem by taking averages over 
the momenta and leaving the remaining dynamics of 
the radiation field hidden in a closure relation fl5ll . In 
addition, systematical approximations are employed to 



reduce the computational cost, either by assuming sym- 
metries and thus reducing the dimensionality, or by av- 
eraging quantities like the energy. 

Practical implementations of transport in compact ob- 
jects have been mostly limited to Newtonian gravity and 
flows jH, EE BH, H3. El HI [H . The motion of photons 
and neutrinos is a genuinely relativistic problem, so ex- 
isting schemes for Newtonian fluid flow often approx- 
imate the coupling terms to order 0(v/c), an assump- 
tion which will be violated as soon as the matter moves 
relativistically itself: in particular, the important case of 
outflows in gamma-ray bursts is excluded in this man- 
ner. General relativistic spacetime curvature can also be 
a substantial effect when compact objects are involved. 

Attempts to discretize the general relativistic trans- 
port problem were mostly focused o n sp herically sym- 
metric spacetimes IMSIMIlElIEIHIlEIITllIllIl, 
since the computational complexity of the problem sim- 
plifies considerably in this case and makes a direct finite 
differencing integration of the Lindquist equation [34] a 
reasonable proposition. 

Discretization schemes not restricted to spherically 
symmetric spacetimes have not been investigated by 
many authors. Schinder j35ll describes the problem in 
the particular case of polar-sliced spacetimes (which at 
that time had been successfully used to simulate the col- 
lapse of a polytrope to a black hole B36ID . Anile and 
Romano 13711 discuss covariant flux-limited diffusion in 
general spacetimes (see also 13^1 ). Cardall and Mez- 
zacappa ll39fl provide a conservative formulation of the 
transport equations, both in the general case and in re- 
striction to spherical symmetry. Farris et al [40] have 
recently added an effective treatment of the radiation 
field to their general relativistic magnetohydrodynam- 
ics code, which discretizes the radiation in the same way 
as the fluid, under the gray approximation and the as- 
sumption that the radiation energy-momentum tensor 
is nearly isotropic in the comoving frame. 

The methods listed above are ultimately based on 
finite-differencing and finite-volume techniques, which 
apply polynomial approximants to reconstruct the radi- 
ation field or moments over it. These schemes have been 
very successful in numerical implementations for hy- 



2 



drodynamics, magnetohydrodynamics, and Einstein's 
field equations. There are at least two major differences 
compared to the radiation transport problem: the first 
is that both the fluid as well as the spacetime metric are 
only dependent on three spatial and one time variable, 
whereas the radiation field depends, at each spacetime 
point, on three momentum variables as well, as men- 
tioned above. The second difference is that the radiation 
field can have strongly peaked distributions in momen- 
tum space which are not easy to approximate with poly- 
nomials. 

This paper investigates a discretization of radiation 
transport based on Monte Carlo methods jilll . Monte 
Carlo schemes are particularly well-suited for high- 
dimensional problems as outlined in section III Bl and 
provide a degree of flexibility which is difficult to 
achieve with polynomial techniques. The scheme is con- 
structed to solve the full transport problem, in general 
spacetimes and including full angular and energy de- 
pendency, but it should be straightforward to restrict to 
particular spacetime symmetries, or to the 0(v/c) ap- 
proximation needed for a coupling to Newtonian flu- 
ids 1 . The scheme is indirect, i.e. it admits an arbitrary 
specification of all probability density functions, which 
is particularly useful to derive adaptive methods for ap- 
plication problems. It can also handle optically thick 
regions by virtue of its hybrid approach, that is by us- 
ing cell-based surface flux estimators. Since the scheme 
is local, it admit parallelization with standard domain- 
decomposition methods. 

The paper is organized as follows: In section [TTJ gen- 
eral relativistic kinetic theory as needed for purposes of 
the Monte Carlo discretization is reviewed, and the ba- 
sic idea of (indirect or deterministic) Monte Carlo tech- 
niques is presented. Section [HI] contains the deriva- 
tion of the Monte Carlo estimators for all terms in the 
collision integral, and also estimators for the energy- 
momentum contained in the radiation field. Section ITVl 
gives an outline of how such a scheme may be real- 
ized and coupled to existing evolution codes for gen- 
eral relativistic magnetohydrodynamics, and discusses 
some particular features of the Monte Carlo approach. 
Finally, section M concludes the paper with a summary 
and short discussion of future work. 



II. BACKGROUND 

A. General relativistic kinetic theory 

Any attempt to represent a gas of photons, neutri- 
nos, or other particles invites the application of kinetic 



theory. The original formalism for general relativis- 
tic kinetic theory is described in the seminal paper by 
Lindquist [34]. This section will review some essential 
aspects of this theory needed for the discrete representa- 
tion, following closely Ehlers |50h , and also Cardall and 
Mezzacappa 13911 . 

In classical physics, the state of an ensemble of n par- 
ticles can be described by a point in the product phase 
space M = Mi x . . . x M n , where each phase space M; 
is spanned by the configuration variables x = z\) 
and the momenta p = (pf/pf, pf)- However, if one is 
rather interested in average properties of the ensemble, 
it is useful to introduce, for each type of particles, a dis- 
tribution function defined on a one-particle phase space: 

Here, N denotes the particle number in a particular 
phase space volume element, i.e. the object / is the 
phase space density of particles of that type. 

A similar object can also be introduced in general 
relativity. Let (M,g) be the space-time manifold with 
metric g. Given the collection T X (M), x G M of tan- 
gent spaces associated with events in M, a new ob- 
ject, the eight-dimensional tangent bundle T(M), can 
be constructed with the metric g ® g and volume form 
- det(g fl6 ) dx° A ... A dx 3 A dp A ... A dp 3 . 

Since particles of mass m and momentum p G T X (M) 
satisfy the mass-shell constraint g{p,p) = —m 2 , the 
(one-particle) phase space of that particle type is given 
byM ffl = {x G M,p e T x (M) r g(p,p) = -m 2 }. 2 If the 
mass shell is parametrized by the three-components of 
the momentum, a volume form on the shell is given by 

n m = dp 1 A dp 2 A dp 3 . (2) 

I Pol 

where po is now a function given by the mass-shell con- 
straint. 

It will later be useful to employ proper frame coor- 
dinates on the mass shell, which we will denote by a 

wedge, i.e. p = (p l ,p 2 ,p 3 ) In this frame the ele- 
ment 7r„, assumes the form 

Tt m = -j^ry dp 1 A dp 2 A dp 3 (3) 



where E(p) = \Jp 2 — m 2 . 

A volume form on the one-particle phase space M m 
can then be constructed canonically, given the configu- 
ration space four-volume form r\ = y ' —g dx° A ... A dx 3 , 



For Monte Carlo schemes as app lied to Newtonian background 
flows, see (HIH@ISIMIiEI3lil and references therein. 



2 This space is sometimes also called the sphere bundle in analogy to 
Riemannian geometry. 
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by the exterior product 



O = Tj A 7T„ 



(4) 



I Pol 



dx° A ... A dx 3 A dp 1 A dp 2 A dp 3 
dx° A ... A dx 3 A dp 1 A dp 2 A dp 3 , 



HP) 

where the second equality is obtained from the specific 
form ©. 

Particle flow in phase space can be approximated by 
separating two different physical regimes: The first are 
long-range interactions with large-scale fields, most no- 
tably the interaction of charged particles with an elec- 
tromagnetic field, and the second are short-ranged in- 
teractions, which we will approximate as point inter- 
actions for the present purpose, in particular contain- 
ing direct interactions with matter (emission, absorp- 
tion, scattering). The former can be represented as a 
four-force, while the latter are modifications to the dis- 
tribution function in any domain containing the interac- 
tion event. 

The particle flow itself is represented by orbits of the 
Liouville operator, which, in the absence of long-range 
forces, is given by the geodesies of the particles and their 
associated momenta. In explicit form, given a connec- 
tion with coordinate coefficients T a jj C , the Liouville op- 
erator is the vector field 



L = 



(5) 



We will often need to integrate over hypersurfaces 
in the manifold M m when considering the particle flow 
through a particular surface. Consider a hypersurface 
E C M m . On E, we can introduce the 6-form [50] 



to 



(6) 



where a a is the vectorial 3-surface element of configura- 
tion space, i.e. o~ a = \ n abcd A dx c A dx d , and n m as 
given above. 

The general relativistic distribution function / is then 
related to the number N[E] of particles crossing the sur- 
face E by the integral 



NIL] 



fcv. 



(7) 



In the specific case where E is a small element of the 
phase space containing an observer at x E M, it can 
be written as E = G x K, where G is a spacelike hy- 
persurface containing x, and K is its associated set of 
mass shells. In the proper frame of this observer, the 
expression then limits (for vanishing spatial volume) 
in the differential f d 3 xd 3 p. Therefore, the function / as 
defined above corresponds to the classical distribution 
function defined in 0. 

From the surface integral 0, one can get a statement 
of particle conservation by considering a domain D of 



the sphere bundle. Application of the surface integral 
to the boundary dD of that domain, and using Stake's 
theorem, one obtains |5(J 



(8) 



AN [3D] = / fco = [ L[f]n, 

JdD JD 



which relates the number of collisions in D to the Liou- 
ville flow L [/] . The differential form of this equation is 
the relativistic Boltzmann or Lindquist equation [34!] : 



l\f] 



dx a 



dp a 



(9) 



and here C [f] stands for a source term which is the limit 
of AN [3D] for a differential domain size D. 

A concept closely related to the distribution function 
is the invariant intensity fl5ll . Starting from the tradi- 
tional definition of the radiation intensity I v , which is 
the energy emitted per time interval, surface area, angle 
and frequency interval, one can construct the Lorentz 
invariant 



(10) 



also called the invariant intensity. This object is related 
to the distribution function by only a numerical factor in 
the Planck constant h and the speed of light c: 



(11) 



The source term C[f] of the Lindquist equation is 
commonly splitted into the invariant emissivity e and 
the invariant opacity a, which in general are functionals 
of / and may contain contributions from spontaneous 
emission, absorption, and scattering: 



C[f]=e[f]-a[f]f. 



B. Monte Carlo techniques 



(12) 



The simplest application of Monte Carlo techniques 3 
is to calculate properties of statistical systems by gen- 
erating numbers with a particular random distribution, 
and thus constructing a sample of the ensemble space. 
This approach is called a direct simulation of a proba- 
bilistic problem, where the probability density functions 
(PDFs) are fixed by the problem to solve. 

Monte Carlo methods can also be applied to determin- 
istic problems, however. In this case, an estimator is de- 
fined, which is a map from the set of samples to a real 
number designed to approximate the value to solve for. 



3 A general introduction to Monte Carlo methods is given in Ham- 
mersley and Handscomb lilll . 
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A simple example is given by the integration problem: 
Consider the definite integral of a function / over the 
real interval [a, b] : 



J a 



dx 



(13) 



A standard numerical technique is to use polynomial 
approximants, e.g. the trapezoidal rule, to obtain a se- 
quence of approximate values which converge to I pro- 
vided / satisfies certain conditions. 

An indirect Monte Carlo method solves this problem 
in a different way: First, random samples X; 6 [a, b], i = 
1, . . . N are generated according to an arbitrary PDF p, 
also defined on [a, b] . Then, the estimator 



{1} 



N 



/fa) 



(14) 



is computed from the samples. For increasing number 
of samples, this estimator approximates the integral I, 
limiting in the correct value for N — > oo, which is also 
expressed by saying that the estimator is unbiased. The 
asymptotic error behaves like 1 / vff. 

A direct comparison to a polynomial technique of 
order k (k = 2 in the case of the trapezoidal rule) 
seems, at first sight, rather unfavorable for the Monte 
Carlo approach: the error of the polynomial scheme 
goes like 1 / N k and is therefore asymptotically superior 
to the Monte Carlo estimator. However, an interesting 
observation concerns multi-dimensional integrals: the 
Monte Carlo estimator retains the error 1 / y/N, whereas 
a polynomial approximation by dimension has an er- 
ror 1 / ~N k / d , where d is the number of dimensions. For 
d = 4, the Monte Carlo scheme and the trapezoidal rule 
have the same accuracy, whereas for d > 4 Monte Carlo 
is asymptotically superior. Therefore, Monte Carlo 
techniques are considered particularly useful for high- 
dimensional problems. 

This fact is of some interest for radiation transport: a 
second-order polynomial scheme should be asymptot- 
ically superior to a Monte Carlo scheme only for four 
dimensions or less. The full radiation transport prob- 
lem, however, is six-dimensional, suggesting that using 
a Monte Carlo scheme may be preferable in that case. 
These estimates are only based on the asymptotic error, 
of course. 

The actual magnitude of the local error can be 
reduced, sometimes by orders of magnitude, with 
variance-reduction techniques Hill like control variates, 
importance sampling and stratified sampling. Im- 
portance sampling, which makes use of the fact that 
the PDF can be chosen arbitrarily in indirect Monte 
Carlo schemes, could be particularly useful for radiation 
transport problems (see discussion in section HV) . 



III. MONTE CARLO APPROXIMANTS FOR 
RELATIVISTIC KINETIC THEORY 

A. Cell discretization, and decomposition of the collision 
integral 

Assume D C M m to be a domain of the phase space 
associated with particles of mass m. From ((8) we see that 
the number of interactions in D is given by the boundary 
integral 



AN[8D] = f fto. 

JdD 



(15) 



Consider configuration space to be parametrized by 
coordinates (t,a,b,c), and a spacetime cell given by a 

region G = [V , V + At'] x . . . x \c\, c\ + Acj], where i,j are 
integers obtained from the discretization. We would like 
to characterize the emission of particles newly created in 
this spacetime interval by means of the source functions, 
the invariant emissivity e and the invariant opacity a. 

For this purpose, first decompose the integral ( fTBI l into 
terms for each surface in configuration space and its en- 
tire associated future-directed mass shells: 



ANldD] 



f CO 



T+ 



fw + ... 



(16) 



where we label the cell surface t — V — const by T , the 
surface t — V + AV = const by T + , and all remaining six 
surfaces likewise. Next, further decompose each surface 
into an ingoing and outgoing component with respect to 
the space-time cell, and label the resulting surfaces Aj , 
A~ ut and so on 4 . 

AN[dD}= f fcv+f fta+( foi+ [ _fw + ... 

JT JT+ J A in J A out 

(17) 

Every point on an outflow boundary surface can be 
backpropagated by the Liouville operator to a point on 
an ingoing boundary surface, and the evolution of the 
distribution function is given by the Lindquist equation 



h 4 



/z 4 df_ 
~? dX 



dj_ 
dX 



(18) 



where A is an affine parameter, and the last relation is 
the usual form involving the invariant emissivity and 
opacity. In general, both e and a are not only functions of 
the local fluid variables (and further macroscopic fields), 
but also functionals of the full distribution function on 
the past light cone at that point, making the relation fl8l 
an integrodifferential equation. 

The numerical scheme presented here separates be- 
tween particles newly emitted in a volume element and 



4 Since all mass shells are future directed, this decomposition is trivial 
for the surfaces T~ and T+. 
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those injected and absorbed subsequently. To formal- 
ize this separation, define, for a particular orbit of the 
Liouville operator, the function ^ to be a solution of 
the pure absorption component with initial condition 
^(Ag) = ^ (Ag) at the ingoing boundary point Ag: 



B. Monte Carlo estimators for emission 



dX 



This function represents the absorption of injected par- 
ticles. Next, define another function J? as the solution 
of the full Boltzmann equation, but for vanishing inflow, 
i.e. ^(Ag) = 0. This represents particles newly created 
(and possibly immediately reabsorbed) inside the cell: 



(20) 



Now construct the sum Jf = J? + J" , which has the 
same inflow condition as ^ , and insert this into the 
Lindquist equation: 



(21) 



-^-=e[f]-a[j?]jf -a[J] 



In the particular case where e and a are independent of 
the changes in the distribution function while propagat- 
ing inside the cell, the function J 1 is a solution to the 
Lindquist equation. 5 

With an expression of the collision integral and the 
separation of ^ into ^ and J> , we define /, /, and / 
by multiplying the corresponding invariant intensities 
by c 2 /h 4 , and decompose (Y7\ into contributions from / 
and /, yielding 

AN[dD] = AN em [dD] + AN abs [dD] (22) 
AN em [dD] : 

AN abs [dD] 



fto + / fto + ... 

T- JT+ 



fto + / fco + ... 

T- JT+ 



Physically, these represent the changes in the particle 
numbers due to emission of new particles as measured 
on the outflow surface, or absorption of particles in- 
jected at the inflow surface. In the next section, we will 
employ a Monte Carlo scheme to solve for the terms 
involving /. The remaining terms are treated in sec- 
tion |Inc] 



5 Note that both the emissivity and the opacity are still allowed to 
depend on the boundary values of ^ and on the affine parame- 
ter A, which in particular admits to model scattering processes. To 
treat the general case in some approximation, it is also conceivable 
to either neglect or approximate the non-linear coupling between 
the partial distributions, and treat each one independently as a non- 
linear ODE. 



All integrals over the ingoing surfaces in (TS3D vanish 
due to the initial condition / — employed in defining 
/. We retain the sum 



(19) AN em [dD] 



f to + / fto 

T+ JA-, 



A out 



fco + ... (23) 



for the effective emission of newly created particles as 
measured on the outflow surface. 

We now define an unbiased Monte Carlo estimator for 
<E3by 

{AN em [dD)} = {[ '?„}+{[ fcv} + {[ fcv} + . 

(24) 

where the estimator for each constituent integral on the 
surface E = T + , A^ uf , ... is obtained from the estimator 

In this expression, N(L) samples have been obtained 
from the evaluation of a probability density function 
P(x,p) on that surface, and the sample values are de- 
noted by (xj, p^. The basic setup of the emission esti- 
mator is illustrated in figured] 

To numerically evaluate the estimators i25) , we intro- 
duce canonical coordinates on each surface (i.e. (a, b, c) 
on T + , (t, b, c) on A~ ut and Ag Ut , and so on) and evaluate 
the form to as in $6$ and |3} 



CO\j+ 



CO\ T - 
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(V A + 



0) 



,b 0Ut 



to\ r - 

*~out 



to\ r + 

'-out 



HP) 

p 

W) 



m 

W) 

HP) 
p h 

HP) 
f 

HP) 
— p c 

HP) 



-gda Adb Adc A dp 1 A dp 2 A dp 3 
-gda Adb Adc A dp 1 A dp A dp 
^gdt Adb Adc A dp 1 A dp 1 A dp 3 
^gdt Adb Adc A dp 1 A dp 1 A dp 3 
-gdt Ada Adc A dp 1 A dp 2 A dp 3 
-gdt Ada Adc A dp 1 A dp 2 A dp 3 
^gdt Ada Adb A dp 1 A dp 2 A dp 3 
^gdt Ada Adb A dp 1 A dp 2 A dp 3 



While the PDF is arbitrary, the factorized form 
P(x, p) = Q(x) R(p) will be useful in practice. The func- 
tion Q(x) can be obtained from a product of uniform 
density functions 

Q(x\x 2 ,x 3 ) = Q^x 1 ) Q 2 (x 2 ) Q 3 (x 3 ) (26) 
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Trapped surface Estimator: 




Figure 1: Emission estimators for an indirect, cell-based Monte Carlo discretization of general relativistic kinetic theory. For 
purposes of illustrating general relativistic effects, the situation near the horizon of a black hole is shown as represented in a 
system of horizon-penetrating coordinates, which are often used in numerical evolutions. Each cell of the discretization scheme 
is bounded by [t, t + At] x [ay a t + Aa] x [bj, bi + Ab] x [c k , c k + Ac], of which only two dimensions are shown here. The light cone 
structure of the spacetime is illustrated on some cell surfaces, and the location of a trapped surface is indicated as well. The orbits 
of the Liouville operator, which are geodesies in the absence of long-range interactions, are indicated as violet curves. Each cell 
surface can be extended to a hypersurface in the tangent bundle M m over the spacetime M, and is then labeled T~, A~ and so 
on. Those subsets which are relevant for the emission estimator are the outgoing subsets of these surface, A~ t and so on, and an 
estimator for the emission part of the collision integral is given by sampling each surface (some samples are indicated by circles) 
and solving the Lindquist equation on the intersection of the orbit with the cell. In general, some of these estimators can also be 
empty, as shown (in this gauge) for AI ut inside the event horizon. 



where the components of x depend on the surface as 
above, and each direction is sampled according to 

Qlix 1 ) = . \ . , (27) 

X high X low 

with x\ , and x\ denoting the boundaries of inteera- 

ntgtl low o c 

tion. 

Having fixed Q(x), we need to specify R(p) as well. 
The use of proper frame coordinates in the momentum 
space volume form n is motivated by the often simpler 
form of emission and absorption coefficients in a frame 
comoving with a fluid. This suggest to use a factor- 
ized PDF R(p) = N(v) M(D.), where v is the frequency 
and Q = (0,0) denotes directional variables in mo- 
mentum space. Sampling frequency and direction sep- 
arately is often useful on physical grounds, and a sim- 
ple form of the directional component of the PDF can be 
obtained from a uniform distribution over the variables 
(cos(0),O). 

There is a complication with selecting suitable sam- 
ples with the form of R(p) given above: Since the in- 
tegrals only contain outgoing points in the phase space 
with respect to that surface, not all samples provide 
valid surface points. It is possible to represent the 
boundary by a transformation from coordinate space to 
a rest frame and adjust the PDF accordingly. However, 
another (simpler) way exists: sample the full directional 
PDF without regard to the boundary surface, but, after 
a transformation of the constructed momentum vector 
to the coordinate frame, reject those samples which are 



ingoing. 6 

Finally, we need a value for / at each sample point. 
This is obtained from the solution of the Lindquist equa- 
tion as in (|20l , under the inflow boundary condition 
^(Ao) = 0. In those cases where e and a are approx- 
imately constants within the cell, we can construct the 
formal solution of the Lindquist equation llT5ll 

/(A) = - + C exp(-aA) (28) 

where C is a constant of integration. W.l.o.g. assume 
Aq = and A = AA at the outgoing point. The condi- 
tion J 1 (0) = leads to C = — ~ and the solution at the 
outgoing point: 

/(AA) = - (1 - exp(-aAA)) (29) 

This is the well-known result that emission in optically 
thick sources (flAA 3> 1) rapidly approaches the source 
function | . If e or a are varying substantially with A in- 
side the cell, equation j20T l can be numerically treated 
with an appropriate ODE integrator. 

Given the estimators ll25l l, we can consider each sum- 
mand a bundle contributing a particular number of par- 
ticles to the total estimate: 

{ k fa,} - m 5 p(*,-,y,-) (30) 



6 This technique, rejection sampling, is generally useful to sample un- 
der constraints, although it can also increase the variance in adverse 
conditions. 
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N(E) 

= EW 

i=l 

These bundles provide a Monte Carlo discretization of 
the distribution function. The estimate {N;} associated 
with a bundle will be called its weight. Physically, we 
may imagine a set of directed bundles of {N,} photons 
to be emitted from the cell. 

In contrast to a sampling of the emissivity following 
a particular fluid world line, which gives a direct Monte 
Carlo scheme and is a grid-less technique as far as the 
radiation field is concerned, the estimator j24l is grid- 
based and indirect, i.e. the PDFs are not fixed and can 
be chosen arbitrarily. We will discuss some particular 
uses of this level of freedom in section HVl 

It is sometimes useful to obtain a representation for / 
itself from the bundles. This can be obtained by assum- 
ing a distributional form 

N(E) 

{/}= £WM(i-ii)^-ft) (3i) 

where x is the tupel of coordinates spanning the partic- 
ular surface, and a, is yet unspecified. We perform the 
integral over the surface E: 



/ 



CO = CO 



N(E) 



EWMCx-XiW-ft) (32) 



i-l 



Here, to has the particular form appropriate for that sur- 
face as above. This allows us to fix the coefficients a; by 
comparison with (f3T~b to 



(33) 



tO(Xi, pi) 

and we therefore obtain 



N(E) 



{/} = E 



.S(x-Xi)S(fi-pi). (34) 



C. Monte Carlo estimators for absorption 

Return to the expression d25l for the emission estima- 
tors: 



co}^ 



n(e; 



N <£ f( x i,Pi) w(*vPi 

1=1 



P(xi,yi) 



(35) 



The particular surface E considered here also appears 
as boundary surface of an adjacent cell D', but with 
a reversed normal vector and thus an inverse associa- 
tion with ingoing and outgoing mass shell components. 
Consider £ to be the complement of E in the union of 
ingoing and outgoing hypersurfaces associated with a 
specific configuration space surface: then both are also 



subsets of dD', and in particular E is ingoing and E is 
outgoing with respect to D'. 

It follows that the estimator l |2"5t also serves as an es- 
timator for the integral over the function / as defined in 
the previous section, effectively giving inflow boundary 
condition in form of a Monte Carlo estimator. A distri- 
butional representation of the function / on E C dD' can 
be reconstructed as in equation l|34l . 

To obtain the propagation and absorption estimator, 
we turn to the second half of the terms in the collision 
integral fTTl i, namely those associated with f: 

AN abs [dD'\ = f fto+ [ fto+ f fco (36) 

Jt- Jt+ Jat„ 



f CO 



Here, all integrals will contribute in the general case. 

Using the distributional representation of / l(3l]l on 
each ingoing surface E,„, we can represent those inte- 
grals by a sum of integrals over non-intersecting, arbi- 
trarily small control regions containing each discretiza- 
tion point (see also figure 

N(E,„) r , 

E 7 ~ \ S(x-x l ) (37) 
■S(p-f)co 

N(E,„) 



{Ni} 

f-[ h, C0(Xi,pi) 



S(x — x l ) S(p — p l ) co 



where the surfaces E; satisfy E; C E !n , E; n Ey = for 
i ^ j, and (x ir pi) £ E,-. 

Now define flux tubes D, from the orbits of the Liou- 
ville operator intersecting E;. These tubes intersect the 
cell at a surface E' contained in the union of its outgoing 
surfaces. Since the tube has been constructed from the 
Liouville flow, the general relation <(8j holds for the tube: 



AN[dDj 



fco 



L[f] n. 



(38) 



Since the tube boundaries are parallel to the Liouville 
orbits, they have zero volume with respect to to = L ■ Cl 
ll50ll and do not contribute to the surface integral. There- 
fore an estimator for the difference in particle number 
between E^ and E; is given by 



(39) 



{AN,} = {if co} - {J J to} 



Introduce a 6+1 decomposition of the tube with the 
generating function given by an affine parameter A. 
Then we can employ the same distributional representa- 
tion as in d3Tt for each slice. The integral over the slices 
defines a function AT; (A) which takes the role of /(A) as 



8 




l ) 6(p — p 1 ) ui 



Figure 2: Absorption estimators: The left panel shows a set of adjacent spacetime cells with intersecting Liouville orbits as in 
figureQ] The cell under consideration is D', for which cells Dj and D2 (among others) provide inflow conditions in the form of 
emission estimators. To obtain collision estimators for the partial distribution f (see text), each sample is extended into a flux tube, 
a seven-dimensional subset of D' generated by the Liouville flow and delimited by small surfaces E,- and E^, as illustrated in the 
right panel. After solving the Lindquist equation for each tube in an averaged sense, the estimators for the outflow components 
of the integrals over / can be constructed from the set of surfaces E^. 



averaged over the slice, and thus follows the same equa- 
tion as / in equation fl9)l : 



dNj(A) 
dA 



-flNi(A), 



(40) 



with initial condition given by {N,}. If a is approxi- 
mately constant inside the cell, we can easily construct a 
formal solution providing the difference estimate: 



{N't} = {Ni} exp(-aAA) 



(41) 



Here AA is the affine distance between the sample 
(xj,pj) and its Liouville development at Hi:. This rela- 
tion reflects the expected result that the particle number 
in a bundle stays constant without source, and decreases 
exponentially with the effective optical depth. In those 
cases where a varies substantially with A, a numerical 
solution for the ODE can be obtained. 

Having obtained an estimate for the number of par- 
ticles on all E' v we can define the remaining part of the 
estimator for / by using equation d3Tb for all outgoing 
surfaces of D: 



{/ /"} = 



N(T, ut) /M'l 



S(x — x 1 ) 



■5{p-p')w 

f-[ Je; 0}{xi,fi 
■6(p-p')u> 



The sample locations are the Liouville developments of 
the samples on the ingoing surfaces, or solutions of the 
geodesic equation in the absence of long-range forces. 
This completes the definition of estimators for all terms 
in the decomposition of the collision integral. 



D. Estimators for energy-momentum 

To couple a radiation transport scheme to Einstein's 
field equations or general relativistic magnetohydrody- 
namics, we need an estimate of the energy-momentum 
tensor R of the radiation field. Since the Monte Carlo 
scheme defines the kinetic function / only in a distri- 
butional sense, the energy-momentum tensor can be es- 
timated in an integral form, which also translates well 
into finite-volume schemes usually applied to obtain a 
correct weak solution for fluids I5U 15211 . 

The starting point for obtaining conservative schemes 
is the general conservation law derived from 33 = 
UK], given a spacetime region Q: 



d*R= / *R = 0. 



(43) 



Given a tangent space basis {e a } and the surface form 
<j a , the surface integral can be expressed in components 
as 



[ *R= I e a R ah 



(44) 



The surface 3H can be decomposed into parts G ; , (J G; = 
3Q as needed for the spacetime discretization. 

The energy-momentum tensor associated with a dis- 
tribution function / is given by [50] 



R 



ab 



K x 



V a P b frt, 



(45) 



where the integration covers the mass shell at the con- 
figuration space point under consideration. 

Inserting this form into equation (|44l , and using the 
distributional representation of / in equation d34b , yields 
the desired estimator for a configuration space surface 

G, c an, 



?a V a V h {/} <?b A 7i 



(46) 
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= J E e a p a {f}u 

M f-[ u>(xi,pi) 
■S(p-Pi)w 

AT(Ei) 

where we have used to = p a a a A n, and where £; is the 
union of G; with all its future-directed mass shells as in 
section iHIBl 



IV. SOME CONSIDERATIONS FOR NUMERICAL 
IMPLEMENTATIONS OF THE SCHEME 

This section contains a brief discussion of how to real- 
ize the Monte Carlo estimators discussed in section HIT1 
in a practical application, and also contains some re- 
marks on potential uses of the arbitrary probability den- 
sity functions in astrophysical models. 

Numerical codes which solve the (magneto-) hydro- 
dynamics equations coupled to Einstein's field equa- 
tions 7 regularly employ a finite-volume discretization 
of the fluid HI HI and either the BSSN 1@ IH HI or 
generalized harmonic [57j| form of Einstein's field equa- 
tions, discretized with finite differences. 

If a radiation transport module is to be coupled to 
such a system, the amount of energy-momentum in the 
radiation field R needs to be accounted for in the BSSN 
or generalized harmonic system, and since conservation 
is satisfied by the total energy-momentum tensor (fluid, 
magnetic field, and all radiation species), all emission, 
absorption and scattering processes extract or supply 
energy and momentum to the fluid. 

A simple approach is to perform the coupling in an 
operator-split manner, in which the radiation field es- 
timators are updated separately. The spacetime projec- 
tion of the cell S employed in section [III] can be identi- 
cal to the cells employed in the finite-volume method, 
but that is not a requirement, and different grid resolu- 
tions can be prescribed as needed. The spacetime met- 
ric updated by the BSSN/harmonic code enters in the 
volume forms and the geodesic equations, whereas the 
fluid variables provide the rest frame and thermody- 
namic state of the matter, possibly including chemical 
and nuclear composition variables, from which emis- 
sion, absorption and scattering coefficients can be de- 
rived. In turn, the estimators <f47b yield approximations 
to the amount of energy-momentum in the radiation 
field and act as sources and sinks for both the spacetime 
and the fluid. 



The radiation transport step itself may be performed 
in this order: 

• Choose a set of probability density functions. 

For each cell surface S, construct a PDF Pe(x, p) 
on the spacetime surface and all associated mass 
shells. This function can be chosen arbitrarily, and 
it can be different between all cells. 

• Create random samples. For each cell surface 
Yj, sample the PDF chosen above with an arbi- 
trary number of samples N(2L). The sample lo- 
cations (xj, on the surface can be obtained by 
generation of pseudorandom numbers in the unit 
interval [0,1], e.g. with the popular Mersenne 
twister algorithm |58] 8 . For a one-dimensional 
PDF, the sample location is then formally the in- 
verse of the cumulative distribution function ap- 
plied to the value on the unit interval obtained 
above llilll . Factorizations like those suggested in 
section UnB] greatly simplify the multidimensional 
case. In general, tabulated values can be employed 
for complicated PDFs. 

• Create new bundles. The samples admit to de- 
fine surface emission estimators as given by equa- 
tion i25h The necessary values of the surface el- 
ement can be obtained from the spacetime metric 
and equations 1261 1. The value of / is in general 
calculated by integrating an ODE, but for constant 
emissivity and opacity a formal solution as dis- 
cussed in section IHI Bl can be used. The effective 
emissivity and opacity are, in the simplest case, 
a function of the material thermodynamic state 
alone, but may also depend on the radiation field 
estimate at inflow when scattering is dominant. 

The affine length AA needed in this expression is a 
function of the intersection of the particle geodesic 
curve with the cell boundaries, which can be ap- 
proximated by a Taylor expansion. 

• Propagate existing bundles. The set of existing 
(and newly created) bundles is propagated up to 
each surface T + . The bundle weight is adjusted ac- 
cording to equation (HQ) , which needs to be treated 
as an ODE in general. In the case of constant opac- 
ity, the affine length can again be approximated by 
a Taylor expansion. 

• Estimate radiation energy-momentum. For 

newly emitted particles, the estimators for the 
energy-momentum in the radiation field, equa- 
tion @7) / approximate the amount of energy and 
momentum removed from the fluid. In turn, 



7 see e.g. |ilHSHHBHBIi,|T3IU0miH 31111 references therein - 



The implementations of the C standard library function rand, or 
its Fortran equivalent, are often not sufficiently accurate for Monte 
Carlo applications. 
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when absorption is estimated, they determine the 
amount returned to the fluid. Scattering is for- 
mally treated in the same way since the end result 
of scattering is always a particular surface state of 
the radiation field, and an associated flux coupled 
to the fluid flux. 

A scheme constructed along these lines covers the 
simulation domain with a large set of bundles and sep- 
arately integrates each bundle. The quality of the fluid- 
radiation coupling, and therefore also the coupling be- 
tween curves in the Liouville flow, is determined by the 
cell size and the number of samples per cell. 

Since the probability density functions are not related 
to the emissivity and opacity, a number of schemes can 
be employed to reduce the discretization error of the 
scheme. Importance sampling jilll is a basic technique in 
which approximately known properties of the solution 
are used to define probability density functions which 
reduce overall variance. In practice, a frequency sam- 
pling roughly based on the local emissivity and opacity 
should give better results than a linear or logarithmic 
distribution function. 

More importantly, the directional distribution of bun- 
dles can be varied freely .This can be used to reduce 
variance in the radiation field in directions which are 
considered more important for a particular problem. 
One example would be the polar heating by neutrino- 
antineutrino annihilation in collapsar disks: If a geo- 
metrical region of interest is identified, the angular PDF 
can be selected such that the field has increased reso- 
lution in those Liouville orbits which intersect that re- 
gion. In essence, the free choice of PDF gives rise to flex- 
ible adaptive techniques, which can be combined with 
actual adaptive mesh refinement in the cell-based dis- 
cretization if required. 

All these advantages can be used in the context of 
Newtonian approximations to the simulations of com- 
pact objects as well: In this case, a suitable background 
metric needs to be derived from the Newtonian grav- 
itational potential, and the coordinate transformations 
to the rest frame are performed to 0(v/c) as in other 
approaches 1151 l23l 15911 . The geometric curvature intro- 
duced by non-Cartesian (spherical, cylindrical) grids is 
already contained in the covariant form of the expres- 
sions. Therefore the scheme is not restricted to general 
relativistic applications. 

The discretization presented here is also conservative, 
which is a desirable property for radiation transport (see 
also 139J). All energy and momentum removed from a 
cell via the flux estimators is accounted for in the expres- 
sion for the energy-momentum tensor, and only absorp- 
tion and out-scattering can reduce the particle number, 
energy and momentum associated with each bundle. 

A final but important remark concerns the paralleliza- 
tion of the scheme. Any full six-dimensional discretiza- 
tion of radiation transport will likely dominate the com- 
putational requirements of a simulation. In recent years, 
one of the traditional approaches to increasing com- 



putational performance, by raising the clock frequency 
of the processing core, has failed due to heat dissipa- 
tion issues l6(ill . Therefore, future architectures will 
tend to employ massive parallelization in the order of 
10 5 . . . 10 6 computing cores or more per machine, and 
numerical algorithms need to be designed to effectively 
use such resources. The hybrid Monte Carlo method is 
particle /cell-based and local, and therefore admits the 
use of standard domain decomposition techniques for 
load distribution. It should be possible to design im- 
plementations which are scalable to large parallel high 
performance architectures in the future. 

V. CONCLUSION AND OUTLOOK 

This paper contains an indirect Monte Carlo dis- 
cretization of general relativistic radiation transport. 
Starting from the general expression for the collision in- 
tegral applied to a cell decomposition of the configura- 
tion and phase space, Monte Carlo functions are intro- 
duced which provide unbiased estimates of phase space 
surface integrals over the relativistic distribution func- 
tion. This approach is indirect, that is it leaves all prob- 
ability density functions used to select samples in phase 
space unspecified. 

The relativistic distribution function is split into two 
components with different inflow boundary conditions 
per each cell. In the linear approximation, where the 
collision terms do not depend on the distribution, this 
decomposition is exact. For practical implementations, 
it provides an approximation which improves with di- 
minishing finite volume of the cell. The split allows to 
treat two parts of the collision integral separately: one 
for the construction of new bundles, and one for the 
propagation of existing ones. In addition to collision es- 
timators, a similar estimator for the energy-momentum 
in the radiation field is constructed, and can be used for 
coupling the radiation to general relativistic magneto- 
hydrodynamics and Einstein's field equations. 

Finally, an outline of a numerical implementation 
of this transport scheme is given, with some addi- 
tional remarks on expectations concerning the adaptive- 
ness, parallel efficiency and conservative nature of the 
scheme. 

The scheme exhibits a number of promising features 
compared to other approaches: the local nature of the 
operators which is expected to help parallelizing the 
scheme, the free choice of probability density functions 
which makes the use of adaptive techniques easier, and 
its fully covariant formulation which should make cou- 
pling to both general relativistic and Newtonian codes, 
possibly with non-Cartesian coordinates, straightfor- 
ward. Besides these properties, an actual implementa- 
tion and the associated tests are expected to be a major 
effort, particularly when a full coupling to general rela- 
tivistic magnetohydrodynamics is considered. 

Some advanced methods may also be worthy of in- 
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vestigation in the future: there exist a number of so- 
called variance-reduction methods, which are ultimately 
designed to reduce the discretization error of Monte 
Carlo estimators. Importance sampling has already 
been mentioned. Quasi-Monte Carlo is the name of a 
set of techniques where sequences are constructed with 
methods from number theory such that they minimize 
a particular distributional measure. Although these se- 
quences are not random, they provide unbiased esti- 
mators with substantially lower variance than random 
sampling. Stratified sampling is another technique which 
cuts the sampling domain into intervals and therefore 
into separate estimators. And control variates split ana- 
lytically tractable control functions from the estimator, 
thus reducing the average discretization error based on 



the approximation made by the variate. More advanced 
techniques are known in the literature, e.g. antithetic 
variables which reportedly reduce variances up to 10~ 6 
in certain cases, according to Hammersley and Hand- 
scomb jilll . 

In the future, Monte Carlo techniques should deserve 
further attention, particularly when full six-dimensional 
radiation transport is to be calculated on upcoming 
large-scale parallel machines. With good parallel scal- 
ing and sufficiently large machine allocations, the simu- 
lation of 10 11 . . . 10 12 bundles should be within the range 
of possibilities during the next five years, which may 
give enough resolution for an interesting class of prob- 
lems. 
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